Introduction aux modèles spatiaux avec R-INLA
François Rousseu
3 mai 2018
Types de données
Il y a 3 principaux types de données spatiales
- Données géostatistiques
- Données de points (point pattern)
- Données d’aires (areal data or lattice data)
Données géostatistiques
Ce sont des données mesurées à des points fixes dans l’espace
- Aussi connu sous le nom de géostatistiques ou kriging.
Données récoltées à des stations météo (température, précipitations, etc.).
geoR, gstat, nlme, glmmTMB, R-INLA
Patrons de points
Ce sont des données où c’est la localisation des points dans l’espace qui est aléatoire
L’ensemble des mentions d’occurrence pour la Phragmite commune (Phragmites autralis).
Un exemple de cas
Nous allons nous concentrer aujourd’hui sur jeu de données concernant le nombre d’espèces de plantes (richesse) à différentes parcelles sur l’île Las Palmas dans les îles Canaries. Cette île est caractérisée par de fortes variations en altitude en raison de la présence de deux volcans sur l’île. Ces données proviennent d’une étude par Irl et al. (2015) et peuvent être téléchargées librement à cette adresse:
https://doi.org/10.5061/dryad.7sk2g/1
Les exemples qui suivent sont fortement inspirés de Zuur et al. (2017).
La base de données légèrement modifiée peut être téléchargée directement à partir de la page GitHub de l’atelier.
d<-read.csv("https://raw.githubusercontent.com/frousseu/introRINLA/master/canary_richness.csv")
head(d)## Plotname easting northing richness altitude TCI
## 1 AC001 221.1035 3190.789 21 0.7784 1.26
## 2 AC002 221.4035 3192.189 19 0.4473 1.14
## 3 AC003 221.8035 3192.689 25 0.2921 1.28
## 4 AC005 222.1035 3192.889 21 0.0610 1.48
## 5 AC006 222.0035 3192.689 23 0.1424 1.78
## 6 AC007 222.4035 3192.989 22 0.1211 1.86On a donc un identifiant de parcelle, des coordonnées, la richesse en espèces et deux variables. L’altitude est en kilomètre et TCI est un index de complexité topographique.
Inspection des données
En premier, utilisons les fonctions du package sp pour convertir le data.frame en objet spatial. Ensuite, nous téléchargeons l’île en question avec la fonction getData du package raster pour la visualisation.
library(sp)
library(raster)
library(scales)
coordinates(d)<-~easting+northing # convert to spatial
proj4string(d)<-"+proj=utm +zone=28 +ellps=WGS84 +datum=WGS84 +units=km +no_defs" # assign projection
ds<-spTransform(d,CRS("+init=epsg:4326")) # transform to latlon
spa<-getData("GADM",country="ESP",level=2) # download Spain as a shapefile
laspalmas<-disaggregate(spa)[221,] # subset Las Palmas
par(mar=c(3,3,0,0)) # plot locations ans Las Palmas
plot(laspalmas,axes=TRUE)
plot(ds,add=TRUE,col=alpha("darkgreen",0.5),pch=16) # alpha adds tranparency to points
On peut également visualiser les localisation avec une carte interactive à l’aide du package mapview.
library(mapview)
mapviewOptions(basemaps = c("Esri.WorldImagery","Esri.WorldShadedRelief"))
mapview(ds,zcol="richness")Avant de faire un modèle, on visualise brièvement les données. Le @data dans le code suivant permet d’extraire le data.frame (la table d’attributs) de l’objet d qui est maintenant un objet spatial, plus précisément un SpatialPointsDataFrame.
plot(d@data[,c("richness","altitude","TCI")])
Un simple GLM
Utilisons premièrement un simple GLM avec un effet quadratique sur l’altitude.
fit<-glm(richness~altitude+I(altitude^2)+TCI,data=d@data,family=poisson)
summary(fit)##
## Call:
## glm(formula = richness ~ altitude + I(altitude^2) + TCI, family = poisson,
## data = d@data)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -4.8273 -1.3622 -0.2198 0.8298 6.8234
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.89824 0.04299 44.159 <2e-16 ***
## altitude 0.58252 0.06330 9.203 <2e-16 ***
## I(altitude^2) -0.58296 0.03322 -17.550 <2e-16 ***
## TCI 0.67431 0.02908 23.186 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 4493.8 on 889 degrees of freedom
## Residual deviance: 2813.5 on 886 degrees of freedom
## AIC: 6655.6
##
## Number of Fisher Scoring iterations: 5On peut visualiser le modèle avec le package visreg. L’argument scale permet d’obtenir les prédictions sous l’échelle de la réponse (et non sous l’échelle log).
library(visreg)
par(mfrow=c(1,2),mar=c(4,4,0,1))
visreg(fit,"altitude",scale="response")
visreg(fit,"TCI",scale="response")
Inspections des résidus
Maintenant, plusieurs des parcelles sont très près les unes des autres et on peut se demander si les nombres d’espèces observés sont plus similaires pour les parcelles près les unes des autres. En d’autres mots, reste-t-il de la variation inexpliquée par nos variables explicatives qui pourrait être structurée spatialement?
Pour ce faire, on peut premièrement extraire les résidus (variation non-expliquée) de notre modèle et les cartographier afin de déterminer visuellement s’il y a des patrons suggérant la présence de corrélation spatiale. Si les petits résidus ou les grands résidus ont tendance à être ensemble, cela nous suggère qu’il y a potentiellement de la corrélation spatiale.
par(mar=c(3,3,0,-0))
plot(d,cex=rescale(resid(fit),to=c(0.2,3)),col=gray(0.5,0.5),pch=16,axes=TRUE)
Variogrammes
Il semble y avoir un tel patron, mais cela n’est pas si facile à confirmer visuellement. On peut vérifier cela de façon plus formelle à l’aide d’un variogramme. Intuitivement, un variogramme représente la variance des différences entre toutes les paires d’observations pour différentes classes de distances.
Pour déterminer ceci, on peut faire un variogramme à l’aide des fonctions du package geoR.
library(geoR)
v<-variog(coords=cbind(d$easting,d$northing),data=resid(fit))## variog: computing omnidirectional variogramplot(v, main="Empirical Variogram for Species Richness",type="b",xlab="Distance (km)",ylab="Semivariance") 
Si les observations près dans l’espace se ressemblent davantage, on devrait s’attendre à ce que les valeurs de variance soient plus faibles pour les distances courtes. En fait, une courbe plate suggère qu’il n’y a pas de corrélation spatiale, alors qu’une courbe qui augmente (et qui atteint possiblement un plateau) suggère qu’il y a corrélation.
Nous avons pris les valeurs par défaut de la fonction variog et notamment le nombre de classes de distances qui est relativement faible. Augmentons le nombre de classes afin d’avoir un peu plus de précision sur ce qui se passe à faible distance.
v<-variog(coords=cbind(d$easting,d$northing),data=resid(fit),breaks=seq(0,20,by=0.5),max.dist=20)## variog: computing omnidirectional variogramplot(v, main="Empirical Variogram for Species Richness",type="b",xlab="Distance (km)",ylab="Semivariance") 
On peut voir que la variance est beaucoup moins élevée à de faibles distances. On peut contrôler davantage le variogramme pour inspecter la forme de la relation.
v<-variog(coords=cbind(d$easting,d$northing),data=resid(fit),breaks=seq(0,10,by=0.5),max.dist=10)## variog: computing omnidirectional variogramplot(v, main="Empirical Variogram for Species Richness",type="b",xlab="Distance (km)",ylab="Semivariance") 
Il semble donc y avoir une importante corrélation jusqu’à environ 2km, après quoi les points sont plutôt plats, indiquant qu’au delà de cette distance, la corrélation est faible ou inexistante.
On a donc une corrélation spatiale et il faut donc en tenir compte dans notre analyse! Pour cela, il faut utiliser des méthodes spéciales qui sont disponibles dans les packages geoR et gstat. Cependant, ces packages peuvent être utilisés avec des modèles gaussiens, mais pas avec les GLM (à part le package geoRglm).
Empiriques vs. Théoriques
Un variogramme empirique est un variogramme généré à partir de données.
## variog: computing omnidirectional variogram
Les variogrammes théoriques sont les modèles qui sont ajustés aux variogrammes empiriques.
library(gstat)
show.vgms(models=c("Exp","Gau","Sph","Cir","Mat"),as.groups=TRUE)
Différents types de variogrammes:
Exp: exponentielGau: gaussienSph: sphériqueCir: circulaireMat: Matérn
Les différents types de variogrammes sont décrits par des paramètres qui peuvent être ajustés.
show.vgms(kappa.range=c(0.5,1,2,5),models="Mat",nugget=c(0.1,0.2,0.3,0.4),max=10,as.groups=TRUE)
Un variogramme est souvent décrit par:
- range (étendue de la corrélation spatiale)
- sill (le plateau atteint en variance)
- nugget (la variation due à d’autres facteurs)
v<-show.vgms(range=2,nugget=3,sill=6,models="Exp",max=10,as.groups=TRUE,ylim=c(0,11),plot=FALSE)
plot(semivariance~distance,data=v[-1,],type="l",ylim=c(0,10),col="blue",xaxs="i",yaxs="i",xlab="distance",ylab="semivariance")
abline(6,0,lty=3)
abline(9,0,lty=3)
abline(3,0,lty=3)
abline(v=2,lty=3)
Modèle spatial
\[y(s) \sim Poisson(\mu(s))\] \[log(\mu(s))=\beta x + \mu(s)\] \[\mu \sim GF(0,\Sigma)\]
Où:
\(y\):
INLA
INLA (Integrated Nested Laplace Approximation) est un algorithme complexe permettant d’approximer des distributions postérieures dans le cadre d’analyses bayésiennes.
Cette méthode peut être appliquée à plusieurs types de modèles que nous utilisons fréquemment (GLM, GLMM, GAM, analyse de survie, etc.). Plus généralement, cette méthode peut être appliquée pour une classe spéciale de modèles, soit les modèles gaussiens latents (Latent Gaussian Models, LGM). Cette classe de modèles généralise plusieurs modèles que nous utilisons couramment.
Pour une introduction à cette approche, je vous suggère Gomez-Rubio (2019). Pour un condensé “accessible” de la théorie derrière INLA, je vous suggère cet article de blogue.
GLM avec INLA
library(INLA)
m<-inla(richness~altitude+I(altitude^2)+TCI,data=d@data,family="poisson")
summary(m)##
## Call:
## c("inla(formula = richness ~ altitude + I(altitude^2) + TCI, family = \"poisson\", ", " data = d@data)")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 0.5112 0.4459 0.1582 1.1153
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## (Intercept) 1.8980 0.0430 1.8138 1.8980 1.9826 1.8978 0
## altitude 0.5828 0.0633 0.4589 0.5826 0.7074 0.5823 0
## I(altitude^2) -0.5831 0.0332 -0.6487 -0.5830 -0.5182 -0.5827 0
## TCI 0.6744 0.0291 0.6169 0.6745 0.7311 0.6747 0
##
## The model has no random effects
##
## The model has no hyperparameters
##
## Expected number of effective parameters(std dev): 4.078(0.00)
## Number of equivalent replicates : 218.25
##
## Marginal log-Likelihood: -3348.80m$summary.fixed[,c(1:3,5)]## mean sd 0.025quant 0.975quant
## (Intercept) 1.8980388 0.04300921 1.8137909 1.9825871
## altitude 0.5827786 0.06333174 0.4588786 0.7074154
## I(altitude^2) -0.5830944 0.03323408 -0.6486835 -0.5182087
## TCI 0.6743690 0.02910015 0.6168778 0.7311290cbind(summary(fit)$coef[,1:2],confint(fit))## Estimate Std. Error 2.5 % 97.5 %
## (Intercept) 1.8982431 0.04298664 1.8142114 1.9827229
## altitude 0.5825224 0.06329847 0.4589618 0.7070942
## I(altitude^2) -0.5829604 0.03321702 -0.6484524 -0.5182393
## TCI 0.6743125 0.02908278 0.6169035 0.7309101Approche bayésienne
Lorsque appliquée aux modèles, on a:
\[P(\theta|data) = \frac{P(data|\theta)P(\theta)}{P(data)}\]
\(data\): données \(\theta\): paramètres du modèle
Priors
Les analyses bayésienne nécessitent des priors. INLA utilise des valeurs par défaut pour les priors associés aux paramètres standards (coefficient d’un GLM par exemple). Ces priors sont définis par une distribution normale centrée sur 0 et ayant une précision (\(\tau\)) de 0.001. La précision est l’inverse de la variance et donc une précision de 0.001 correspond à un écart-type de:
\[ \sigma = \sqrt{\frac{1}{\tau}} = \sqrt{\frac{1}{0.001}} \approx 31.6\]
En général, il faut s’assurer que ces priors sont appropriés pour notre analyse. Pour l’intercept, la précision utilisée est de 0. Les détails sont expliqués dans ?control.fixed.
SPDE
SPDE veut dire Stochastic Partial Differential Equation. Cette approche se base entre autres sur une discrétisation du Gaussian Random Field (GF ou GRF) par un Gaussian Markov Random Field (GMRF). La fonction utilisée pour décrire la covariance des observations est la fonction de Matérn.
GLM spatial
Étapes
L’estimation d’un modèle avec un effet spatial avec INLA requiert de passer par plusieurs étapes.
- Mesh: permet de créer une triangulation qui discrétise le GF
- Projector Matrix: permet de relier les localisations à la grille
- SPDE: permet détablir le lien le SPDE et la fonction Matérn
- Spatial Field: permet de créer le champ spatial
- Stack: permet d’intégrer les différents éléments
- Formula: formulation du modèle
- Modèle: permet de faire tourner le modèle
Voici un résumé graphique de ces différentes étapes avec les fonction associées.
Mesh
La première étape consiste à créer une grille (mesh) qui va être utilisée pour approximer le champ Gaussien. En premier, il faut obtenir les localisations.
locs<-coordinates(d)
head(locs)## easting northing
## 1 221.1035 3190.789
## 2 221.4035 3192.189
## 3 221.8035 3192.689
## 4 222.1035 3192.889
## 5 222.0035 3192.689
## 6 222.4035 3192.989Par la suite, on utilise la fonction inla.mesh.2d pour créer la grille. Optionnellement, on peut utiliser une grille non-convexe qui respecte davantage le contour de nos observations.
library(INLA)
hull<-inla.nonconvex.hull(locs,convex=-0.05)
mesh<-inla.mesh.2d(loc=locs,offset=c(1,5),max.edge=c(1,3),cutoff=1,boundary=hull)par(mar=c(0,0,6,0))
plot(mesh,asp=1)
points(locs,col=alpha("black",0.6),pch=16,cex=0.7)
offset
Spécifie l’étendue de l’extension de la grille au delà des observations et l’extension d’une zone tampon au-delà de cette grille. Cette dernière permet de réduire les effets de bordure lors des estimations. L’étendue de cette zone tampon doit être au moins aussi grande que l’étendue de la corrélation spatiale.
max.edge
Permet de spécifier les dimensions maximales des triangles de la grille à l’intérieur et à l’extérieur dans la zone tampon. Plus ces triangles sont petits, plus les approximations sont précises, mais plus les calculs sont longs. En général, il n’est pas nécessaire que les triangles dans la zone tampon soient aussi petits qu’à l’intérieur de la grille.
cutoff
Par défaut, chaque observation sera utilisée comme un coin d’un triangle (vertex). Pour éviter la création de trop nombreux petits triangles, on spécifie une valeur à cutoff au delà de laquelle des points voisins seront ignorés lors de la création des triangles.
boundary
Cet argument permet d’utilier l’étendue non-convexe autour des observations.
Il faut éviter d’avoir des triangles avec des aigles trop aigus. Autrement, les estimations sont de moins bonnes qualités.
Projector Matrix
Ceci permet de relier les localisations à la grille et d’établir une pondération qui permet d’estimer les valeurs du champ spatial pour chaque localisation. Les localisations situées sur les coins (vertex) seront estimées à partir des valeurs de la grille et les localisations à l’intérieur du triangle seront estimées à partir d’une moyenne pondérée calculée en utilisant les trois coins du triangle dans lequel elles se trouvent.
A<-inla.spde.make.A(mesh,locs)Dans cet exemple, l’observation en rouge est située à l’intérieur d’un tiangle et les valeurs représentent la pondération qui sera utilisée sur les valeurs de chaque coin pour estimer la valeur de ce point.
par(mar=c(0,0,5,0))
plot(mesh,asp=1,xlim=c(223,226),ylim=c(3172,3176))
points(locs,pch=16,cex=2)
i<-889 # ligne contenant la localisation en rouge
points(locs[i,,drop=FALSE],col="red",pch=16,cex=2)
w<-which(A[i,]>0) # pondérations associées au point
points(mesh$loc[w,1:2],label=1:nrow(mesh$loc),pch=1,cex=4,lwd=2,col="red")
text(mesh$loc[w,1:2],label=round(A[i,w],2),adj=c(-0.6,0.5),font=2,col="red")
SPDE
Ceci permet de définir les éléments du SPDE et les éléments associés aux caractéristiques du champ spatial.
spde<-inla.spde2.pcmatern(mesh,alpha=2,prior.range=c(2,0.5),prior.sigma=c(5,0.01))L’argument alpha est pour spécifier une des paramètres de la fonction de Matérn. Ce paramètre doit être fixé entre 0 et 2. Nous devons également définir les priors du champ spatial. L’approche par défaut est d’utiliser des quantiles pour définir les priors sur l’étendue (range) et l’écart-type (standard deviation, sd) associés au champ spatial. Cette façon de faire se base sur l’approche des Penalized-complexity priors développée par Simpson et al. (2017). En résumé, cette approche permet de pénaliser l’étendue de la corrélation spatiale vers l’infini (ce qui réduit la complexité du phénomène) et la variance vers 0 (ce qui réduit également la complexité du phénomène). En pratique, le prior pour l’étendue \(r\) de la corrélation est spécifié avec \(\alpha\) et \(r_{0}\):
\[P(r<r_{0})= \alpha\]
Où \(\alpha\) représente la probabilité que \(r\) soit inférieur à \(r_{0}\). Dans notre cas, prior.range=c(2,0.5) veut dire que nous croyons qu’il y a 50% de chances que le l’étendue de la corrélation spatiale soit supérieure à 2km et donc qu’il y a également 50% des chances q’elle soit inférieure à 2km. Le prior pour la variance du champ spatial, \(\sigma\), est spécifié avec \(\alpha\) et \(\sigma_{0}\):
\[P(\sigma>\sigma_{0})= \alpha\]
Où \(\alpha\) représente la probabilité que \(\sigma\) soit inférieur à \(\sigma_{0}\). Dans notre cas, cela indique qu’on croit qu’il y a 1% des chances que la variance spatiale en richesse soit supérieure à 5. Ici, il faut nous rappeler que nous travaillons sous l’échelle log et que \(e^5\approx148\) espèces.
Spatial Field
Ceci permet de définir un index qui sera utile pour spécifier ou récupérer les éléments associés au champ spatial.
spatial.index<-inla.spde.make.index(name="spatial",n.spde=spde$n.spde)Stack
Le stack est une façon compliquée de fournir les données, les variables et les effets à INLA. Ceci n’est pas essentiel pour des modèles simples, mais ça le devient lorsque les modèles sont plus compliqués ou lorsque l’on veut par exemple générer des prédictions à partir de notre modèle.
d$altitude2<-d$altitude^2
v<-c("TCI","altitude","altitude2")
X<-data.frame(Intercept=1,d@data[,v])
stack<-inla.stack(data=list(richness=d$richness),A=list(A,1),effects=list(spatial=spatial.index,as.list(X)),tag="est")Formula
On créé la formule décrivant le modèle souhaité. Notez que l’intercept est créé à la mitaine afin de.
model<-richness~-1+Intercept+altitude+altitude2+TCI+f(spatial,model=spde)Modèle
Finalement, on fait tourner le modèle en appelant la fonction inla et en fournissant les différents arguments. La spécification de compute=TRUE fait en sorte que les distributions postérieures seront calculées pour toutes les observations. La spécification de config=TRUE permettra plus loin de faire des simulations et de générer des valeurs aléatoires à partir de notre modèle.
m<-inla(model,data=inla.stack.data(stack),control.predictor=list(A=inla.stack.A(stack)),family="poisson")
summary(m)##
## Call:
## c("inla(formula = model, family = \"poisson\", data = inla.stack.data(stack), ", " control.predictor = list(A = inla.stack.A(stack)))")
##
## Time used:
## Pre-processing Running inla Post-processing Total
## 1.0284 7.4206 0.4006 8.8496
##
## Fixed effects:
## mean sd 0.025quant 0.5quant 0.975quant mode kld
## Intercept 1.9971 0.1018 1.7965 1.9973 2.1967 1.9976 0
## altitude 0.6313 0.1694 0.2989 0.6311 0.9646 0.6307 0
## altitude2 -0.5258 0.0804 -0.6837 -0.5259 -0.3678 -0.5259 0
## TCI 0.4624 0.0511 0.3618 0.4625 0.5623 0.4628 0
##
## Random effects:
## Name Model
## spatial SPDE2 model
##
## Model hyperparameters:
## mean sd 0.025quant 0.5quant 0.975quant mode
## Range for spatial 2.1249 0.3606 1.4863 2.1026 2.9009 2.0639
## Stdev for spatial 0.4277 0.0272 0.3772 0.4267 0.4841 0.4245
##
## Expected number of effective parameters(std dev): 216.05(7.575)
## Number of equivalent replicates : 4.119
##
## Marginal log-Likelihood: -2894.36Visualisations
Le modèle généré (m) est un objet complexe contenant énormément d’information où il n’est pas toujours facile de s’y retrouver. Heureusement, comme nous venons de le voir la fonction summary peut être utilisée pour un sommaire rapide. Si on veut extraire plus d’informations sur notre modèle, la première étape est probablement d’inspecter les différents éléments de la liste formant m.
names(m)## [1] "names.fixed" "summary.fixed"
## [3] "marginals.fixed" "summary.lincomb"
## [5] "marginals.lincomb" "size.lincomb"
## [7] "summary.lincomb.derived" "marginals.lincomb.derived"
## [9] "size.lincomb.derived" "mlik"
## [11] "cpo" "po"
## [13] "waic" "model.random"
## [15] "summary.random" "marginals.random"
## [17] "size.random" "summary.linear.predictor"
## [19] "marginals.linear.predictor" "summary.fitted.values"
## [21] "marginals.fitted.values" "size.linear.predictor"
## [23] "summary.hyperpar" "marginals.hyperpar"
## [25] "internal.summary.hyperpar" "internal.marginals.hyperpar"
## [27] "offset.linear.predictor" "model.spde2.blc"
## [29] "summary.spde2.blc" "marginals.spde2.blc"
## [31] "size.spde2.blc" "model.spde3.blc"
## [33] "summary.spde3.blc" "marginals.spde3.blc"
## [35] "size.spde3.blc" "logfile"
## [37] "misc" "dic"
## [39] "mode" "neffp"
## [41] "joint.hyper" "nhyper"
## [43] "version" "Q"
## [45] "graph" "ok"
## [47] "cpu.used" "all.hyper"
## [49] ".args" "call"
## [51] "model.matrix"Distributions postérieures
Les premiers éléments qui devraient nous intéresser sont les distributions postérieures des coefficients associés à nos différentes variables. On peut extraire les quantiles de celles-cià partir du sommaire des effets fixes.
m$summary.fixed## mean sd 0.025quant 0.5quant 0.975quant
## Intercept 1.9971341 0.10184491 1.7964590 1.9973040 2.1966707
## altitude 0.6313364 0.16944080 0.2989431 0.6311220 0.9645868
## altitude2 -0.5258156 0.08042319 -0.6837075 -0.5258632 -0.3678136
## TCI 0.4624125 0.05107728 0.3617522 0.4625407 0.5622671
## mode kld
## Intercept 1.9976386 8.520360e-07
## altitude 0.6306963 5.063866e-07
## altitude2 -0.5259492 4.770584e-07
## TCI 0.4628019 2.483898e-06On peut également représenter graphiquement ces distributions à l’aide des distributions marginales complètes. Comme on peut le voir, tous les coefficients sont relativement loin de zéros.
par(mfrow=c(2,2),mar=c(4,4,1,1))
invisible(
lapply(names(m$marginals.fixed),function(i){
p<-m$marginals.fixed[[i]]
plot(p[,1],p[,2],type="l",xlab=i,ylab="density")
abline(v=0,lty=3)
})
)
Paramètres spatiaux
Le sommaire des paramètres associés au champ spatial s’accèdent par le sommaire des “hyperparamètres” (hyperparameters). Avec INLA, les paramètres fixes représentent les paramètres de la régression et tous les paramètres de variances ou associés au champ spatial sont considérés comme des hyperparamètres.
m$summary.hyperpar## mean sd 0.025quant 0.5quant 0.975quant
## Range for spatial 2.1249450 0.36058394 1.4863228 2.1025889 2.9008522
## Stdev for spatial 0.4277459 0.02724472 0.3771544 0.4267242 0.4840661
## mode
## Range for spatial 2.0638723
## Stdev for spatial 0.4245108On peut également illustrer les distributions postérieures de ces paramètres.
par(mfrow=c(1,2),mar=c(4,4,1,1))
invisible(
lapply(names(m$marginals.hyperpar),function(i){
p<-m$marginals.hyperpar[[i]]
plot(p[,1],p[,2],type="l",xlab=i,ylab="density")
})
)
Champ spatial
Le champ spatial peut être visualiser de différentes façons. La façon la plus simple est d’utiliser les fonction de INLA pour projeter les valeurs du champ sur une grille. Ceci peut être fait de cette façon.
library(fields)
library(viridisLite)
# https://haakonbakka.bitbucket.io/btopic108.html#92_plotting_the_spatial_mean_field
xlim<-range(d$easting)
ylim<-range(d$northing)
# - Can project from the mesh onto a 300x300 plotting grid
proj<-inla.mesh.projector(mesh,xlim=xlim,ylim=ylim,dims=c(300,300))
# - Do the projection
mfield<-inla.mesh.project(projector=proj,field=m$summary.random[['spatial']][['mean']])
sdfield<-inla.mesh.project(projector=proj,field=m$summary.random[['spatial']][['sd']])
par(mfrow=c(1,2),mar=c(3,3,2,3))
image.plot(list(x=proj$x,y=proj$y,z=mfield),col=viridis(100),asp=1)
axis(1)
axis(2)
image.plot(list(x=proj$x,y=proj$y,z=sdfield),col=viridis(100),asp=1)
axis(1)
axis(2)
Prédictions
Les prédictions dans INLA sont générées en soumettant des observations pour lesquelles la variable réponse est NA (dans notre cas, richness=NA). Si on utilise l’argument compute=TRUE, INLA se chargera de calculer des distributions postérieures pour l’ensemble des valeurs avec NA. Cela peut nous permettre par exemple de soumettre des valeurs d’altitude pour étudier son effet sur la richesse.
On génère 50 valeurs d’altitude entre les valeurs minimales et maximales.
n<-50
x<-seq(min(d$altitude),max(d$altitude),length.out=n)
newX<-data.frame(Intercept=1,TCI=mean(d$TCI),altitude=x,altitude2=x^2)Pour l’exercice, on va assumer que l’effet spatial est constant en prenant un endroit quelconque sur l’île. On répète cette localité autant de fois que le nombre de valeurs d’altitude. Idéalement, il serait plus logique de ne pas tenir compte de l’effet spatial dans les prédictions, mais cela complexifie passablement le code nécessaire.
newlocs<-matrix(c(224,3175),ncol=2)[rep(1,n),]On associe la localisation fictive à la grille.
A.pred<-inla.spde.make.A(mesh=mesh,loc=newlocs)On construit une stack pour fournir les valeurs à prédires et on la joint à celle générée plus haut qui contient les données en tant que tel. On spécifie tag=pred ce qui va nous permettre d’identifier les lignes contenant les prédictions.
stack.pred<-inla.stack(data=list(richness=NA),A=list(A.pred,1),effects=list(spatial=spatial.index,as.list(newX)),tag="pred")
stack.full<-inla.stack(stack,stack.pred)Avec ce tag=pred, on construit un index qui nous permettra de récupérer les lignes contenant les prédictions.
index.pred<-inla.stack.index(stack.full,tag="pred")$dataFinalement, on fait tourner notre modèle en prenant soin de spécifier compute=TRUE pour que INLA calcule les distributions postérieures. On spécifie également link=1 pour s’assurer que INLA retransforme les valeurs sous l’échelle du nombre d’espèces et non sous l’échelle log.
m<-inla(model,data=inla.stack.data(stack.full),control.predictor=list(A=inla.stack.A(stack.full),compute=TRUE,link=1),family="poisson")On récupère les valeurs prédites à partir du sommaire des valeurs prédites (summary.fitted.values) en utilisant l’index créé plus haut.
p<-m$summary.fitted.values[index.pred,c("0.025quant","mean","0.975quant")]
plot(x,p[,"mean"],ylim=range(unlist(p)),type="l",xlab="Altitude en km",ylab="Richesse")
lines(x,p[,"0.025quant"],lty=3)
lines(x,p[,"0.975quant"],lty=3)
Liens utiles
Livres
Advanced Spatial Modeling with Stochastic Partial Differential Equations Using R and INLA by Krainski et al. (2019). Disponible en ligne.
Bayesian inference with INLA and R-INLA by Gómez-Rubio (2019). Disponible en ligne.
Une liste de livres sur http://www.r-inla.org/books
Articles importants
Ces articles sont généralement très avancés et difficiles à comprendre. Le premier est peu peu plus accessible que les autres et est un bon point de départ. Les livres mentionnés plus haut et les différents exemples sur le site de R-INLA sont généralement de meilleurs sources pour se familiariser avec R-INLA.
Spatial modeling with R‐INLA: A review by Bakka et al. (2018)
Bayesian Computing with INLA: A Review by Rue et al. (2017)
Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations by Rue, Martino et Chopin (2009)
An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach by Lindgren, Rue et Lindström (2011)
Site et forum de R-INLA
Blog
Une introduction “accessible” à INLA (mais tout de même assez théorique!):
A gentle INLA tutorial